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Corrective effect of many-body interactions in dynamical 
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ABSTRACT 

Dynamical friction is a fundamental and important phenomenon in astrophysics. The 
Chandrasckhar formula is a well-known analytical estimation of the effect. However, 
current astrophysicists have realized that the formula is not correct in some cases 
because of several approximations dared in the formulation and/or complex non- 
linearities in the real universe. For example, it has been indicated that the dynamical 
friction doesn't work in cored density profiles (constant density in the central region) 
despite that the Chandrasckhar formula predicts drag force even in the constant densi- 
ties. In the former half of this paper, I discuss by N-body simulations that many-body 
interactions are also important in actual dynamical friction though derivation of the 
Chandrasekhar formula is based on the assumption of two-body interaction. In the 
simulation, the many-body interactions are caused by a very small number of field 
particles co-rotating with a perturber. However, the contribution from the many-body 
interactions accounts for a non-negligible fraction of the actual dynamical friction. 
In the latter half, I discuss why the cored profiles suppress the dynamical friction. 
One possible explanation is that corrective effect of the many-body interactions drive 
orbital motion of the perturber. The cessation of dynamical friction by this corrective 
effect would be feasible even in shallow cusp density profiles although the shallow cusp 
may evolve into a constant density. 
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1 INTRODUCTION 

Dynamical friction (hereafter, DF) is an important physical 
process for a wide variety of dynamics, e.g. stellar dynam- 
ics, galaxy and planetary formation etc. That's why it is in 
great demand to seek deeper u nderstanding of t he na ture 
of DF. As a historical milestone. IChandrasekhatl (|l943l ) has 
succeeded in formulation of the effect of DF. However, sev- 
eral approximations were unavoidable in the derivation of 
the formula. The Chandrasekhar formula was derived un- 
der the assumptions of uniform density and isotropic ve- 
locity field. Moreover, interaction between a perturber and 
field particles was supposed to be integration of countless 
two-body interactions and orbital periodicity of all particles 
was not taken into account. By efforts of many numerical 
studies, current astrophysicists have realized that the Chan- 
drasekhar formula is not perfe ct and loses accuracy of the 
estimation in some cases (e.g. iTremaine fc Weinberg! 1 1984 
IWeinberdH989l : I Jiang fc Binnevll2000l ; iFuiii et al.ll2006t ). 



In this paper, I concentrate on the suppression 
of DF on a globular clus ter (hereafter, GC) in a 



Goerdt et al.l 2006: Sanchez-Salcedo et all 2006 


: Read et al. 


2006 


; lAngus & Diaferid 120091: ICowsik et all 


20091: llnoue 



20091 ). However, the discussion in this paper can be gen- 
eralized to other astronomical bodies by scaling physical 
units. The drag force by DF is negligibly weak for GCs in 
large systems, like our Galaxy. In contrast, in small sys- 
tems like dwarf galaxies, the drag force is far too strong (see 
iBinnv fc Tremamell2008l . chap. 8). Thus, the GCs in dwarfs 
are presumed to lose their orbital energy and fall into the 
galactic centre by strong friction force in time-scale of the or- 
der of ~ 1 Gyr (|Tremainc 1976; He rnandez fc Gilmore|[l998 ' 
Oh et all |2000|; IVesperinil |2000|. l200ll: iGoerdt et all 12006 
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dwa rf galaxies even now and are as old as the age of universe 
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to as 'the dynamical friction problem' 
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However, a solution of the problem has been proposed. 
If dwarf galaxies have a cored dark matter halo which has a 
constant density region in its centre, the DF against the halo 
is weakened considerably, and GCs should be able to survive 
beyond the age of the universe. By an analytical approach 
using the Chandrasekhar formula, iHernandez fc Gilmord 
(1 19981 ) discovered that a King-model halo can weaken the 
DF in the core reg i on. Recently , by N -bod y simu l ations , 
iGoerdt et ail (|2006h . iRead et all (|2006h and llnouel (120091 1 
confirmed complete cessation of DF in a core region of 
haloes. Thus, although the Chandrasekhar formula indicates 
gradual shrinkage of the GC orbit, the N-body simulations 
demonstrate complete defunctionalization of DF. This dis- 
crepancy between analytic and numerical studies (the Chan- 
drasekhar formula and N-body simulations) still remains. 
Besides, the cessation of DF is a topic not restricte d to 
dwarf galaxies and GCs. iGualandris fc Merrittl (|2008l ) has 
discussed the same core-stalling of DF on binary black hole 
against a cored stellar cluster and suggested that fluctuation 
of the halo potential weakens the DF in the core. 

My goals in this paper are to clarify hidden physics 
in the nature of DF and the mechanism of the DF sup- 
pression in cored density profiles by N-body simulations. 
Constant density distribution is called, in other words, 'har- 
monic oscillator potential', which means that all particles 
in the potential have the same orbital period. In this sense, 
harmonic potentials are expected to induce complex effects 
in non-linear manner and be beyond the applicable domain 
of the Chandrasekhar formula. It has been known for a long 
time that DF shows curio u s beh aviors in constant densi- 
ties. ITremaine fc Weinberg! (1 19841 ) analytically argued that 
a constant density core causes energy feedback on a per- 
turber by dynamical resonance in positive or negative sense 
de pending on the density distri b ution. In the sam e context 
as ITremaine fc Weinberel i| 19841 ). IWeinberd l| 19891 ) indicated 
that density response of the density field to a perturber or- 
biting in the core region can exert torqu e on the perturbe r 
and drive the orbital motion. Moreover, IBoilv et al.l (|2008l ) 
also analytically argued that in the harmonic potentials the 
perturber acts as a catalyst for the evolution of the energy 
distribution of the field particles by the orbital resonance. 

By the way, the cored dark haloes are inconsis- 
tent with results of cosmological N-bod y simulations (e.g. 
iNavarro et all 1 19971 : ISpringel et al .1120081 ). However, among 
small scale galaxies, many observations of low surface bright- 
ness (LSB) galaxie s have confirmed the cor e d dark haloes 
(e.g. lBurkertill995l; ISalucci fc Burkertll2000l ; |de Blok et al.l 
l200ll : lLi fe Chenll2009l ; lde Blokll2010l . and refernces therein), 
and nearby dwarf galaxies also favor the cored haloes 
(|Gilmore et al.ll2007l ). In addition, numerical simulations in- 
cluding baryonic physics also predict core creation by bursty 
super novae in the centre of dwarf haloes in cosmologica l con- 
text (|Mashchenko et alj 120081; iGovernato et al.l koiOl ) and 



isolated models ( Pasetto et al.|[2*010l ). 



This paper is organized as follows: I will explain my 
simulation method and present my simulation results in §2. 
In §3, I will show analyses of the data generated by my 
simulations. In §4, I will summarize my results and present 
discussion. 



2 THE SIMULATIONS 

My numerical simulation is a typical pure N-body simula- 
tion ( no gas component). I use Barnes-H ut modified tree- 
code (|Barnes fc HutJ 1 19861 : iBarnesI Il990l ) , setting an open 
angle of 6 — 0.5. A special-purpose calculator for collision- 
less N-bod y simulations, GRAPE-7, is used with the tree 
algorithm ( Makinol Il99ll ). The leapfrog time-integrator is 
adopted with the shared-timestep method. The total num- 
ber of timesteps is 11841 for the whole of a simulated period 
which corresponds to 10 Gyrs in the real time-scale. 



2.1 The settings 

In this paper, I use Burkert profile (|Burkertlll995l ) as a cored 
halo model of dwarf gal axies, which is motivated from obser- 
vatio ns of LSB galaxies (|Salucci fc Burkertll2000l : lLi fc Cher] 
12009^ and numerical simulat ions of formation of dwarf galax- 
ies (|Mashchenko et al"1l2006l ). 



p(r) 



Porl 



(r + r )(r 2 + rg) ' 



(i) 



po is the effective core density, ro is the core radius which 
defines the constant density region. In my simulation, I set 
po = 0.1Mq/pc 3 , ro = 1.0 kpc, and the virial mass of 
the halo M mr = M{r vir ) = 2.0 x 10 9 M Q , the virial ra- 
dius r V ir = 9.88 kpc. Eq[T] estimates the mass enclosed in- 
side 300 pc to be M(300pc) = 8.79 x lO 6 ^ Thi s value 
is consistent with the result of IStrigari et al.l ()2008h which 
has suggested that all dwarfs have the same mass scale, 
Af(300pc) ~ 10 7 Mq. I add an exponentially decaying en- 
velope to prevent instability at the outer region caused by 
the artificial cut-off radius, rw (|Springel fc White]|l999l ). 
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and c = r V i r /ro- I set r decay = O.lrw. With this envelope, 
total halo mass enclosed in r cutoff = 3.16 x 10 4 pc becomes 
Mtotai = 2.44 x 10 9 M o . 

Velocity dispersion is given by the solution of Jeans 
equations, 



a r (r) 



p(r) 



dr'r p(r') 



(3) 



where /3 is the anisotropy parameter. In this paper, I assume 
the isotropic velocity state in the halo, setting — (oy = 
c"e = cr<t>)- The velocity distribution is determined by local 
MaxweUian approximation, 

F ^ = 4 K2^) 3/VeXP (-2^)' ^ 

wher e F(v) is a prob ability distribution function of veloci- 
ties (|Hernouisdli993l ). By the assumptions of isotropic and 
MaxweUian velocity distribution, I satisfy the suppositions 
underl ying the derivatio n of the Chandrasekhar formula, al- 
though iLhTetal] (|2005l ) has reported slight instability in the 
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Figure 1. The time-evolutions of orbital radius of the GC. The Figure 2. Density profiles of field particles of analytic model 

Chandrasekhar formula is calculated with In A = 3.72. (EqQ]andr2]l and in the N-body simulation for t=0 and 10 Gyr. 



Burkert profile realized by the local Maxwellian approxima- 
tion. 

I conduct N-body simulations by 10 7 particles to re- 
solve the halo. All halo particles have the same mass, m = 
2.44 x 10 2 M Q , with a softening length of e = 3.0 pc. Before 
operating actual simulations, the halo model needs to be re- 
laxed to some extent. I run a simulation without GCs for 2 
Gyrs. After the relaxation, I employ the relaxed halo as the 
initial state (t — 0) of the halo. 

In my simulation, the GC is represented by a point 



mass with m„ 



2.0 x 10 Mq. 



the GC particle is set to e g 



The softening length of 
10 pc. Initial orbital ra- 



dius of the GC is set to 750 pc. The initial orbit is cir- 
cular, and I employ a single GC in my simulation. Multi- 
G C cases are al s o inte resti ng, but have been investigated 
by iGoerdt et al.1 (|2006f ) and llnouei (|2009h . which didn't in- 
dicate any large differences from the single GC cases. If the 
GC is resolved by many particles, some additional effects 
may play important roles: mass-loss, dynamical heating, etc 
(e.g.lFuiii et al.ll2006l: iMiocchi et al]l2006l: lEsauivel fc Fuchl 
120071 ; ICapuzzo-Dolcetta fc Miocchil 120081 ). Besides, stellar 



comp onents in a dwarf galaxy can enhanc e the DF of 
GCs (jSanchez-Salcedo et al.|[200g ; lBekkill2010l ). However, in 
this paper, I concentrate on investigations of fundamental 
physics of DF, and do not consider such effects. 



2.2 The results 

The time-evolution of orbital radius of the GC is plotted in 
Figd] As a comparison, I also plot an analytic calculation 
by the Chandrasekhar formula, 



~dt 

dr 
dt 



4-7T In AG 2 pm gc 



erf(X) - ^e- 



r dv c 
v c dt 



(5) 



(6) 



where X = v c /(\/2a). v c is circular velocity of the GC as a 
function of radius. G is the gravitational constant. For the 
Coulomb logarithm, I set In A = 3.72, which is determined 
by the best-fitting to the result of simulation for t < 4 Gyr. 
If I take the original definition, InA = \n{h max /b m in) = 
ln(ro/e) = 5.81, where b max , m i n are the maximum and the 



minimum impact parameter. In the derivation of Eq|6l the 
assumption of circular orbit is used ( Hernandez fc Gilmord 
1 19981 ; Ilnouell2009l ; ICowsik et all 120091 ). As seen from FiglTl 
in the N-body simulation, the orbital shrinkage stops after 
the GC entered into the constant density region (r < 300pc) . 
On the other hand, the Chandrasekhar formula fails to re- 
produce the N-body result in the core region, in spite of the 
excellent fit outside the core. In this paper, I define t < 5 
Gyr as 1 exerting DF phase', and t > 5 Gyr as ''suppressed 
DF phase'. 

The density profiles of the halo are plotted in Fig[2] 
Even in the end-state of the simulation, the density profile 
is stable. The energy conservation rate of the system, 1 — 
E end /E ini , is 1.92 x 10~ 4 . 



3 THE ANALYSIS 

3.1 The exerting DF phase 

In the N-body simulation, for t < 5 Gyr, DF is clearly ex- 
erted on the GC, leading to orbital shrinkage. The Chan- 
drasekhar formula with InA = 3.72 succeeds in estimation 
of the effect in this phase. 

Firstly, I calculate dynamical energies of all particles. 
The specific energy of the ith particle is 



Ei 



2 



Gm„ 



E 



Gm 



+ e 2 



(7) 



where £ij, gc mean position vectors of the i, jth field particles 
and the GC, respectively. Next, I calculate increase of Ei at 
intervals of 1 Gyrs. In Fig(3] I show a distribution of all field 
particles in E-AE diagram. AE corresponds to the increase 
of Ei between 2 and 3 Gyr. In this figure, I can find a notable 
feature. At a certain energy, — log | E/$o\ ~ —0.05, there 
are some particles of which energies changed largely. This 
feature looks like two horns. Thus, I define — log|_E/$o| < 
-5.23 x 10~ 3 and AE/\$ \ > 1-27 x 10~ 2 (inside the red 
square) as positive-horn (P-horn), — log |i5/«I>o| < —5.23 x 
10" 3 and AE/\$ \ < -1.27 x 10" 2 (inside the green square) 
as negative-horn (N-horn). The numbers of P- and N-horn 
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Figure 3. Contour plot of distribution of field particles in a 
E — AE diagram. The ordinate indicates an amount of increase 
of dynamical energy of each particle between 2 and 3 Gyr: 
AE = E(3Gyr) — E(2Gyr). The abscissa indicates energy before 
the change: E(2Gyr). The values are normalized by |$rj| which 
is absolute value of the potential energy at r = ro analytically 
estimated. Colour-bar labels mass in a bin, log(M/Mo). Black 
colour means no mass. In this diagram, I separate the particles 
into three groups; positive- horn (the red square), negative- horn 
(the green square), the others. 
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Figure 4. Positions of the P-, N-horn particles and the GC at 
t = 2 Gyr. I set X-axis to the direction of the GC from the 
galactic centre. Z-axis coincides with the GC angular momentum 
vector. Y-axis lies in the GC orbital plane. The GC is rotating 
counter-clockwise in the X-Y plane. 



particles in this plot are 1,513 (3.69 x 10 5 Mq in total) and 
2,028 (4.95 x 1O 5 M in total), respectively. 

In Fig[J] I show positions of the P- and N-horn particles 
at t=2 Gyr. The particles are projected on the GC orbital 
plane. By calculating dot products of angular momentum 
vectors of the GC and the horn particles, I confirm that 
most of the particles plotted have positive dot products at 
t=2 Gyr: prograde rotations with the GC (1,194 of 1,513 
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Figure 5. Examples of typical orbits of the P- and N-horn parti- 
cle in rotating coordinate with the GC. Although the horn parti- 
cles are defined by energy increase in t=2-3 Gyr, the orbits plotted 
here are for t=l-4 Gyr. The large (small) dots indicate positions 
of the particles at t=l Gyr (t=2 and t=3 Gyr). 



P-horn and 1,734 of 2,028 N-horn particles have prograde 
rotations). Thus, I find that most of the P-horn particles 
are orbiting behind the GC, whereas, the N-horn particles 
are orbiting ahead of the GC. Moreover, since most of these 
particles have short distances from the GC and the prograde 
rotations, I can ascribe the horn-feature in Fig[3]to the simi- 
larity between orbits of the horn particles and the GC. These 
particles strongly interact with the GC. Conversely, the GC 
ought to be affected by these particles. 

I examine orbital motions of the horn particles. In Fig[5] 
I illustrate typical orbits of the P- and N-horn particles 
in the coordinate rotating with the GC. The P-horn par- 
ticle in this figure approaches the GC from behind. The 
orbit expands when catching-up, and finally retreats from 
the GC after the orbital expansion. This orbital expansion 
corresponds to the energy increase seen in Fig(3] On the 
other hand, the orbit of the N-horn particle is the reverse, 
comes close to the GC (in inertia coordinate, the GC is 
catching up to the N-horn particle), shrinks the orbital ra- 
dius, and runs away from the GC. As a result, the energy 
of the N-horn particle decreases. The mechanism of these 
orbital behaviors is simple. Under density gradients, gener- 
ally, inner orbits have shorter orbital periods, on the other 
hand, outer ones have longer periods. The P-horn particle 
is pulled forward by gravity of the GC, and gains energy. 
However, this energy increase causes the orbital expansion 
and makes angular velocity slow. In the case of the N-horn 
particle, the particle is pulled backward and decrease the 
energy. This energy decrease results in the orbital shrink- 
age and faster angular velocity. This phenomenon, that the 
particle slows down in azimuth when pulled forward and 
speeds up when held back, has come to be called 'the don- 
key e ffect' (|Lvnden- Bell fc KaInaiJl972l ; lBmnv fc Tremamel 
120081 ). The orbits in Fig[5]imply that higher energy particles 
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Figure 6. Cumulative energy transfers to the GC from the P- 
and N-horn particles in the exerting DF phase. The short-dashed 
line corresponds to the sum of both. Increase in Et means energy 
injection into the GC, decrease means energy absorption from the 
GC. 



lose orbital energy, on the other hand, lower energy particles 
gain energy. This behavior can be confirmed in Fig[3] which 
indicates that the two horns lean slightly from top-left to 
bot tom-right. 

iFuiii et ah! (|2009h has also investigated this kind of or- 
bital motion (Trojan horseshoe orbits) in N-body systems , 
which look like the orbits in Fig[5] (see IFuiii et all 120091 . 
Fig. 7). It is worthy of special mention that these orbits ex- 
hibit quite different behavior from two-body interaction that 
is assumed in the derivation of the Chandrasekhar formula. 
This orbital behavior of the horn particles is presumed to 
be caused by many-body interaction, rather than two-body 
interaction (see, Appendix A). 

To probe the energy transfer from the horn particles 
to the GC, I rerun the same simulation while calculating 
mechanical work on the GC, 



E T {t) 



Vg C (t') ■ a Bpc (t')dt' . 



(8) 



This integration means cumulation of the integrand on each 
timestep. a spc represents net acceleration exerted on the GC 
by a specific particle group (e.g. the P-, N-horn or the other 
particles) . 

I compute energy transfers on the GC in the same sim- 
ulation and show the result in Fig(6] Though the horn par- 
ticles were defined by energy change during 2-3 Gyr, I draw 
the energy transfers during t=1.5-3.5 Gyr in Fig(6] setting 
tint = 1-5 Gyr in Eq(8] In agreement with a naive expec- 
tation, I find that the P-horn particles absorb energy from 
the GC and the N-horn particles inject energy into the GC. 
It is notable that the behaviors of the injection and absorp- 
tion are monotonic. The transfers are sharp in 2-3 Gyr, slow 
in t < 2 Gyr and t > 3 Gyr. This implies that the energy 
transfers are transient. 

The energy transfers in Fig[5] don't completely bal- 
ance; however, the energy absorption by the P-horn particles 
slightly surpasses the injection by the N-horn. I compare the 
energy transfer on the GC by both of the P- and N-horn par- 
ticles with that by the other particles in Fig[7] The figure 
indicates that the other particles also steal energy from the 
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Figure 7. Comparison of cumulative energy transfers to the GC. 
The solid line labelled as 'horn particles' is the same as 'sum' in 
Fig|6]but setting ti n ; = 2 Gyr in Eq[8] 



GC. This energy loss of the GC seems to be caused by DF 
against the other field particles except the horn particles. 
Interestingly, the energy loss caused by the horn particles is 
comparable to that by the other particles, accounting for a 
half of actual DF by all particles. The dotted line labelled as 
'Chandrasekhar formula' is calculated by substituting Eq[5] 
for a 3 p C in Eq(S] using the data of GC orbit in the simula- 
tion. The DF from all particles is well fitted by the Chan- 
drasekhar formula with In A = 3.72 This result implies that 
actual DF can be separated into two modes. The first mode 
is the well-known conve ntional DF based on the concept 
of IChandrasekharl (|l943l ). This is a DF that is attributed 
to random encounters with a myriad of field particles in an 
isotropic velocity distribution. On the other hand, unlike the 
first mode, the cause of the second mode is a small number 
of particles which are in many-body interaction; the horn 
particles. The horn particles have only 0.0354 per cent of 
the total mass of the system. 

However, keep in mind that the criterion for defining 
the horn particles is somewhat arbitrary. Moreover, the re- 
sult in Fig[7] doesn't necessarily mean that DF caused by 
'true 1 many-body interacting particles accounts for a half of 
the actual DF. Even in two-body interaction, if a field par- 
ticle has a similar velocity vector to Vg C and a small impact 
parameter, the particle causes strong interaction and large 
energy transfer. There would be no clear-cut way to extract 
only the true many-body interacting particles. However, it 
is sure that the horn particles defined in Fig(3] include many 
many-body interacting particles. 



3.2 The suppressed DF phase 

For t > 5 Gyr, the DF doesn't seem to act on the GC to 
shrink the orbit. This phenomenon itself has already been 
known. iHernandez fc Gilmord (|l99fet ) found that this core- 
stalling of DF can be explained by the Chandrasekhar for- 
mula. However, their estimation was based on the assump- 
tion of completely constant density in the core. Without 
this assumption, the Chand rasekhar formula predicts slow 
but certain orbita l shrinkage (|Goerdt et al j2006l ; lRead et alj 
l200fj| ; llnouell2009r i. as seen in Figffl 
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Figure 8. The same figure as Fig|3] but the differential is taken 
between 7 and 8 Gyr: AE = E(8Gyr) - E(7Gyr). Criteria for 
defining P-horn (N-horn) particles is — log \E/$o\ < —3.0 X 10~ 2 
and AE/\$ \ > 1.27 x 1CT 2 (A£/|* | < -1-27 x 1CT 2 ). 
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Figure 9. Positions of P-, N-horn particles and the GC at t = 7 
Gyr. 



On the other hand, the N-body simulations indicate the 
complete suppression of the DF in the cored region. The 
Chandrasekhar fo rmula is no longer viable in this phase. 
iRead et al.l (|2006T ) has suggested that the mechanism of 
the DF suppression is anisotropic velocity distribution in- 
duced by the GC. They attributed the DF suppression to in- 
crease of number of particles in 'co-rotating state', in which 
more field particles rotate with the GC in p rograde direc- 
tion than initial isotropic state. Ilnoud (|2009t ). however, re- 
vealed that the increase of co-rotating particles is marginal 
and anisotropy induced by the GC is very weak. Besides, 
iGualandris fc MerrittJ (|2008l ) suggested that fluctuation of 
the halo potential weakens the DF in the core. But, they 
supposed the case of a black hole and a stellar cored struc- 
ture, the perturber being so heavy that their theory is not 
applicable to the case of a dwarf galaxy and GCs. 

I execute the same analyses as in the previous subsec- 
tion. In Fig[8] I redefine P- and N-horn particles by the 
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Figure 10. Examples of typical orbits of the P- and N-horn parti- 
cle in the coordinate rotating with the GC. Although the particles 
are defined by energy increase in t=7-8 Gyr, the orbits plotted 
here are in t=6-9 Gyr. The large (small) dots indicate positions 
of the particles at t=6 Gyr (t=7 and t=8 Gyr). 



same way. The numbers of the P- and N-horn particles are 
613 and 834, respectively. Positions of these horn particles 
at t=7 Gyr are plotted in Fig[9] As well as the exerting DF 
phase, I confirm that most of the particles plotted have pro- 
grade rotations with the GC (399 of 613 P-horn and 761 of 
834 N-horn particles). 

In Fig llOl I illustrate typical orbits of P- and N-horn 
particle. This figure indicates very complex orbits. Both par- 
ticles come close to the GC from the front, come inside the 
GC orbital radius, rotate around the GC, expand the or- 
bits and finally go away behind the GC. It seems as if the 
particles dodge the GC. Although these particles are respec- 
tively defined as a P- and a N-horn particle, these orbits are 
identical. Regardless of the complexity, by considering ab- 
sence of the donkey effect in the constant density region, 
the following explanation can be made. In constant densi- 
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Figure 11. Cumulative energy transfers to the GC from the P- 
and N-horn particles in the suppressed DF phase (t; n ; = 6.5 Gyr 
in this plot). 



ties, orbital period becomes independent of orbital radius. 
In coming close to the GC, the particles are orbiting slightly 
outside the core with slower angular velocities than the GC. 
By the gravity from the GC, the particles are pulled back- 
ward and lose energy. The particles sink toward the centre 
by the energy loss. However, by the absence of the donkey 
effect in the core region, the angular velocities of the parti- 
cles are kept constant even in the inner radii. After getting 
behind the GC, the particles are pulled forward by the GC 
and gain energy. As a result, the orbits of the particles ex- 
pand again and are exiled from the core. This orbit is also 
not in two-body interaction but in many-body interaction. 
In Appendix A, I will discuss about this orbit in restricted 
three-body problem. 

In Fig llll I show cumulative energy transfers into the 
GC. In consistency with the complexity of the orbits, the 
mechanical works on the GC also show very intricate fea- 
tures. The transfers are no longer monotonic, unlike in the 
exerting DF phase (Fig(6|. The N-horn particles continue 
to inject energy until t ~ 8 Gyr, but absorb energy from 
the GC afterwards. Meanwhile, the P-horn particles inject 
a little energy at first, but largely absorb later. Thus, the 
total transfer from all horn particles indicates an intricate 
waving behavior. Although these particles are defined by 
E(8Gyr) — E(JGyr), the particles largely inject and absorb 
energy in t < 7 Gyr and t > 8 Gyr. These behaviors implies 
not transient but long-term interaction between the particles 
and the GC. 

From the discussion above, in the suppressed DF phase, 
it is clear that the energy transfers from some particles have 
long duration. In this case, I couldn't extract all many-body 
interacting particles successfully. In order to overcome this 
difficulty, I screen particles in t=5-6, 6-7, 7-8, 8-9, 9-10 Gyr 
by calculating AE. If a particle is defined as a horn particle 
even only once in the five time-periods, I regard the particle 
as a many-body interacting particle in all time-periods. I 
gather all many-body interacting particles in this five time- 
periods. The numbers of particles listed up are shown in the 
top of Table [T] While operating the same simulation start- 
ing from t=5 Gyr, I compute energy transfers into the GC 
from the gathered horn particles and the others. The re- 



Table 1. The numbers of particles defined as horn particles in 
each time-intervals. The top and bottom tables are lGyr-interval 
and 0.5Gyr-interval, respectively. The leftmost column indicates 
time-period in which I calculate AE to define the horn particles. 
The bottommost row is the number of particles for which all par- 
ticles listed above are gathered. Because some particles overlap, 
the number of particle gathered is not identical to simple sum of 
rightmost columns. 
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6 - 7 Gyr 
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Figure 12. Cumulative energy transfers to the GC from the gath- 
ered horn and the other particles in the suppressed DF phase by 
setting the time-interval to calculate AE to be 1 Gyr. For com- 
parison, an analytical estimation by the Chandrasekhar formula 
is also plotted, but it is reduced by half from the naive estimation 
by Eq[5] 



suit is shown in Fig |12l In this figure, the cumulative energy 
transfer from the many-body interacting particles gradually 
increases. On the other hand, the other particles absorb en- 
ergy from the GC. These cumulative transfers, however, in- 
dicate periodic flapping. This periodicity appears just in 1 
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Figure 13. The same plot as Fig |12l but using 0.5Gyr-interval 
to extract the gathered horn particles. 



Gyr period, implying an artificial effect caused by setting 
the time-interval to calculate AE to be 1 Gyr. 

Therefore, I follow the same procedure to calculate AE 
with 0.5Gyr-interval. The numbers of particles listed up are 
shown in the bottom of Table [1] I show cumulative energy 
transfers in Fig |13l with an analytical estimation by Chan- 
drasekhar formula for the sake of comparison. However, the 
analytical estimation is halved, according to the result of 
§3.1. The periodic flapping in Fig |12l disappears. Moreover, 
the energy injection from the horn particles are monotonic 
for t < 9 Gyr. The energy absorption by the other particles 
is also monotonic. Interestingly, the behavior of the energy 
transfer from the horn particles is largely different from that 
in the exerting DF phase; the horn, many-body interacting, 
particles continue to inject energy into the GC in the sup- 
pressed DF phase. On the other hand, the other particles 
absorb energy from the GC as well as in the exerting DF 
phase. The halved analytic formula provides excellent fit to 
the absorption for i < 9 Gyr though there is a slight devia- 
tion at t ~ 8 Gyr. These influences from horn and the other 
particles are canceled out and the GC orbit can remain on 
a constant orbital radius as if the DF on the GC vanished. 



4 DISCUSSION AND SUMMARY 

By the N-body simulation in this paper, I conclude as fol- 
lows. In the exerting DF phase, the Chandrasekhar formula 
seems good at estimating the actual DF. Nevertheless, quan- 
titatively, a very small number of particles accounts for a 
non-negligible fraction of the actual DF although they oc- 
cupy less than 0.1 per cent of the halo mass. These parti- 
cles are in many-body interaction, which have quite differ- 
ent orbits from two-body interaction. The presence of such 
m any-body interacting p articles disagrees with the concept 
of IChandrasekharl (|l943l ). However, in his formula, uncer- 
tainty of determining In A would conceal the contribution 
from these particles. 

In the suppressed DF phase, the many-body interacting 
particles alter direction of energy transfer and feed energy 
into the perturber. On the other hand, the other particles 
exert DF, the orthodox DF, ceaselessly. In this phase, as 
a result, these energy transfers cancel out with each other. 



The orbit of the perturber is thus stabilized. I suggest that 
this could be a possible mechanism of the suppressed DF in 
a constant density structure. 

Although my results indicate a failure of the Chan- 
drasekhar formula, this is not surprising. The derivation 
of the formula has assumed an unrealistic circumstance: 
uniform density, isotropic velocity field and non-periodic 
straight orbit. Moreover, a perturber in a real galaxy has to 
rotate around center of the galaxy. Therefore, from the very 
beginning it is unreasonable to apply the Chandrasekhar 
formula to all kinds of density distributions. 



iTremaine fc Weinberg! <|l984h and IWeinberd (|l989h 
have analytically suggested energy feedback on a per- 
turber by densi ty response of field. Especially, figures in 
IWeinbera (|l989h are very intelligible schematically. He in- 
dicated shapes of density wakes. In his results, outside core 
region, the density wake is asymmetric and exerts drag force 
on the perturber. On the other hand, inside the core, the 
wake become symmetric and the drag force disappears. The 
symmetric shape is similar to the distribution of the horn 
particles in my Fig(5] although the horn particles in my sim- 
ulation is selected by energy response, not density. However, 
in my simulation , the symmetric density wake indicated in 
IWeinberd (|l989h can not be observed in the suppressed 
DF p hase. Besides, the N-body simulations of I Read et al.l 
(|2006h . almost the same simulations as mine, also don't in- 
dicate the sy mmetric density w ake, but asymmetric wake 
(see Fig.4 of iRead et al.l ( |2006h ) even inside the core re- 
gion. iGoerdt et al.l (|201C ) also showed similar asymmetric 
desity wake in their Fig. 5. This discrepancy may or may 
not be cau s ed by low resolution in our simulations. Also, 
Boilv et al. I (|2008h has suggested that the field particles on 
circular orbits lose or gain energy according to the relative 
phase to the orbit of the perturber while the perturber plays 
the role of a catalyst to modify the distribution function of 
the field particles and induces anisotropy. 

The result of this paper indicates that not only two- 
body interactions but also many-body interactions are re- 
sponsible for the actual DF. The many-body interacting 
particles, horn particles, are very small in number and co- 
rotating with the perturber. Moreover, the resistance to 
the DF by this corrective effect is also feasible in shallow 
cusp density profiles although the shallow cusp may evolve 
into a constant density (see, Appendix A). This implies 
sensitivity of the actual DF to velocity field of the halo 
models. I employed in this paper the isotropic assumption, 
P — 0.0 in Eq(3j as the most fundamental model. However, 
many numerical simulations have sugges ted dependence of 
f3 on radius in dwarf and other galaxies (Ivan ^lbaddll982l ; 
iMashchenko et al.ll2008l ; lBinnv Tremainell2008h . These re- 
sults imply that actual DF against a live halo would depend 
sensitively on the distribution functions of the haloes. In ad- 
dition, there are some generalised applications of th e Chan- 
drase k har formula ; again st aspherical dist ributions ([Binnev 
1977 ; iPesce et all I1992T I. vari able In A (iHashimoto et al. 
20031 ). DF on spherica l bodies jEsquivel fc Fuchdl2007l) . in 



a gaseous medium(e.g. lOstrikerl 19991 ; iKim et al.ll2008l ). etc. 
The results in this paper alert to thoughtless usages of the 
Chandrasekhar formula. 
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APPENDIX A: CIRCULAR RESTRICTED 
THREE-BODY PROBLEM IN THE CORE 
REGION 

I discuss here why and how the orbits indicated in Fig llOl 
appear in the suppressed DF phase. Since the orbit of the 
perturber is nearly circular in my simulation and density 
field is stable, I expect that circular restricted three-body 
problem would be a valid method to investigate the orbits 
in the core region. Here, I assume that the GC is orbit- 
ing at a separation of -Ro = 300pc from the galactic cen- 
ter in a circular orbit. Let the enclosed mass of the halo, 
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Figure Al. Contour plot of effective potential in the case where 
the perturber is rotating inside the core region with ro = Ikpc. 
Colour- bar labels the energy, & e ff/\&o I- The coordinate center 
is the center of mass of the system. The dot-dashed, solid and 
dashed white lines indicate equipotcntial surfaces of Lagrangian 
points LI, L2 and L3, respectively. The green line indicates an 
isosurface on which the orbits indicated in Fig llOl are possible. 



Mo = 4-7T p(r)r 2 dr. I set the coordinate center to the 
center of mass of the system. Then, the galactic center and 
the GC are located at xm = [— Rom gc /(Mo + m gc ), 0, 0] and 
x gc = [RqMq/(Mo + m gc ), 0,0], respectively. Then, I cal- 



culate effective potential at x 
rotating with the GC, 



[x, y, 0] in the coordinate 



m 

Rgc 



(x +V) 



(Al) 



where Rm = \x — xm\, R gc = \x — x gc \. The first and second 
term are gravitational potential of the halo density distribu- 
tion and the GC, respectively. The third term is centrifugal 
potential. 

In Fig lAll I show a contour plot of the effective poten- 
tial in the case of that the core radius is set to ro = Ikpc, 
which corresponds to the suppressed DF phase in my simula- 
tion. I also plot equipotential surfaces of Lagrangian points, 
LI, L2 and L3. The green line indicates a equipotential sur- 
face between the effective potentials at L2 and L3: $ e //,£2 
and $ e //,£3- This green line runs through around the GC 
just like the horn-particle orbits in Fig llOl If field particles 
have a Jacobi integral higher than <& e //,z,2 and lower than 
$ e //,£3 5 some of these can orbit along this green line. I con- 
sider the particles shown in Fig llOl to be examples of this 
kind of particle. 

Special attention should be paid to the fact that the 
effective potential at L2 (the solid white line) is lower than 
that at L3 (the dashed white line). This height relation be- 
tween <& e //,L2 and $ e //,£3 permits appearance of orbits like 



Figure A2. Contour plot of effective potential in the case of cusy 
density distribution with ro = 50pc. The dot-dashed, solid and 
dashed white lines indicate equipotential surfaces of Lagrangian 
points LI, L2 and L3, respectively. The green line indicates an 
isosurface on which the orbits indicated in Fig[5]are possible. 



Fig llOl On the other hand, in cases of cuspy density distri- 
bution, $ e // l L2 is generally higher than $ e ff,L3 (see below). 

In Fig |A2l I show a contour plot of the effective poten- 
tial in the case where the core radius is set to ro = 50pc, 
which corresponds to the case of cuspy density distribution: 
the exerting DF phase. In this case, $ e //,i2 is higher than 
&eff.L3 - Then, the green lines in this case can not run across 
between potential minima, i.e., no orbit can run through be- 
tween the GC and the galactic center since such a route is 
completely prohibited. The green lines in Fig |A2l correspond 
to the horn-particle orbits in Fig[5] 

In the case of the cored density, if the GC is orbiting 
in even inner region of the core, much more particles can 
become horn-particles shown as the green line in Fig lAll 
Therefore, force against DF becomes strong toward inner 
region, at a certain radius, the force balances with the DF 
on the GC. This means that the DF suppression is not at- 
tributed anisotropic velocity state ind uced by the GC ro - 
tation : co-ro tating state suggested by iRead et al. | (|2006t) . 
Ilnouei (|2009l ) has shown that such a co-rotating anisotropic 
velocity state is very marginal. If the appearance of horn- 
particles is the true mechanism of DF suppresion, even 
though the velocity distribution is isotropic, core-stalling of 
DF is possible. Although the horn-particles I suggest are 
co-rotating particlles in a sense, but there is no need to in- 
duce the anisotropy in the velocity field. This is consistent 
wi th the margi nally anisotropic velocity distribution found 
bv llnouel (|2009t ). 

Next, I address the cases of shallow cusp density distibu- 
tions by means of the circular restricted three-body problem. 
Some studies have concluded that actual dark matter profile 
of dwarf galaxies would have a shallow inner cusp with slopes 
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Figure A3. Sectional effective potential surface on x-axis (on the 
line of the GC and halo centre) in the case of 7 = 0.3. 
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Figure A4. Sectional effective potential surface on x-axis in the 
case of 7 = 0.6. 



between 0.3-0.6 jGovernato et al.ll2010l ; lOh et al.ll20lbT ). In 

order to investigate the cases of shallow cusps, I use another 
density profile, 



p(r) = 



Po 



(r/r )T[l + (r/r y 



|(/3-7)A 



(A2) 



where po = O.IMq /pc 3 and ro = l.Ofcpc. I set the parame- 
ters, a = 1.5, f3 = 3.0 (outer slope) and 7 = 0.3 or 0.6 (inner 
slope). If a = 1.5, /3 = 3.0 and 7 = 0.0, this profile becomes 
very similar to Burkert profile. In this profile, I assume that 
the GC with the mass of m gc = 2.0 x 10 Mq is orbiting 
at a separation of Ro — 300pe from the galactic center in a 
circular orbit. 

As I discussed above, some field particles can be in the 
special orbit (Fia llOp if the effective potential on L3 point is 
higher than that on L2. This fact implies that all I have to 
do is to calculate the effective potential on L3 and L2. Then, 
I plot the effective potential along x-axis (the line of the GC 
and the halo centre). Fig |A3l shows the effective potential in 
the case of 7 = 0.3. In this case, the effective potential on L3 
is higher than that on L2. I expect that some horn particles 
in the special orbit would be existent. Then, according to 
the main discussion of the paper, these horn particles could 
counter and weaken DF on the GC. 

I show the effective potential in the case of 7 = 0.6 in 
Fig |A4l The position of the GC, the scale radius, ro, and 
other parameters are the same. The effective potential on 
L3 is still higher than that on L2. However, the difference 
between these is very small in comparison with the previous 
case. Only a minority of particles can be in the special orbit 
if 7 = 0.6. In this case, the horn particles may be unable to 
balance against the DF on the GC. 

From the discussion above, steeper inner slope of the 
density profile allows fewer particles to orbit in many-body 
interaction manner like Fig llOl However, in this discussion, 
it is very difficult or impossible to investigate how many 
particles in this orbit are required to cancel out the DF 
on the GC. After all, detailed inves tigation needs to o per- 
ate N-body simulations. However, as lRead et al.l (|2006l ) has 
concluded, N-body simulation suffers from the decline of the 
inner cusp. Then, I can not conclude what value of 7 is the 
criterion of the DF stalling in (nearly) cored density profiles. 



Since steeper (shallower) cusp allows fewer (more) horn par- 
ticles, the number of particles required to observe this kind 
of orbits in N-body simulations would depend on the in- 
ner slope. In the cases of steeper cusps, to operate N-body 
simulation needs more particles. 



